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Choptuik has demonstrated that naked singularities can 
arise in gravitational collapse from smooth, asymptotically 
flat initial data, and that such data have codimension one in 
spherical symmetry. Here we show, for perfect fluid matter 
with equation of state p = p/3, by perturbing around spheri- 
cal symmetry, that such data have in fact codimension one in 
the full phase space, at least in a neighborhood of spherically 
symmetric data. 
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1. Cosmic censorship and critical collapse 

Spacetime singularities are a ubiquitous feature of gen- 
eral relativity. In astrophysical situations, singularities 
only arise inside black holes. While no information (light 
signal) can leave the black hole region of spacetime, the 
singularity itself cannot emit any light signal at all: It 
is spacelike, that is, it can be thought of as a local end 
to time. Nevertheless, the field equations of general rel- 
ativity also admit solutions with timelike singularities. 
Such a singularity may be loosely thought of as a point 
in space, just like a particle, and in many solutions light 
signals could travel from such a singularity to infinity. 
The spacetime to the future of the singularity is then no 
longer uniquely determined by initial data in the past. It 
is therefore an interesting question to ask how generic or 
natural such spacetimes, or the initial data from which 
they evolve, are. A "cosmic censorship conjecture" has 
been formulated by various authors, stating roughly that 
naked singularities do not arise from asymptotically flat 
data for reasonable matter The only counterexam- 
ples existing until recently could be written off as either 
not asymptotically flat, involving matter models (dust) 
that form singularities even in the absence of gravity, 
or requiring high symmetry (spherical shells). Until re- 
cently, however, it could still be conjectured that no 
smooth, asymptotically flat data for a reasonable matter 
model would evolve into a naked singularity. That was 
disproved by the spacetimes constructed numerically by 
Choptuik p|, which we review now. 

Consider a smooth one-parameter family of spherically 
symmetric, asymptotically flat, smooth initial data for a 
self-gravitating scalar field. The only condition on the 
family is that it contains both data sets which form a 
black hole (say, for large parameter p), and data sets 
which do not (say, for small p). Choptuik showed that 



by fine-tuning the value of p to a critical value p„, one can 
obtain arbitrarily small black holes. For p P* the black 
hole mass scales approximately as -Mbh — C(p — p*) 7 , 
with 7 ~ 0.37 a universal constant. A solution evolving 
from fine-tuned initial data approaches one and the same 
solution, independently of what the initial data looked 
like. The better the fine-tuning, the longer the solution 
follows this "critical solution" . The critical solution is 
an attractor of codimension one, because it has precisely 
one growing perturbation mode. Its basin of attraction 
is the black hole threshold, because the eventual fate of 
the solution depends on the amplitude of the growing 
mode in the following manner. If it is present with, say, 
positive amplitude, the final outcome will be a black hole. 
If it is present with negative amplitude, the final outcome 
will be dispersion of the scalar waves to infinity, without 
a black hole. If its amplitude is precisely zero, in the 
limit of infinitely good fine-tuning of the parameter p, 
the solution settles down to the critical solution and never 
leaves it. 

The crucial point for cosmic censorship is that the 
critical solution has a naked curvature singularity. Any 
smooth, asymptotically flat initial data in which the 
growing mode is not present approach the critical so- 
lution and its naked singularity. These data form a set 
of codimension one (in spherical symmetry) because only 
the amplitude of the one growing mode needs to be fine- 
tuned to zero. The spacetimes arising from data on the 
black hole threshold do not contain a "zero mass black 
hole" , but contain a region of the universal critical solu- 
tion surrounding the naked singularity, which is smoothly 
matched to an asymptotically flat region that depends on 
the initial data. Near the singularity the deviation of the 
actual collapse spacetime from the universal critical solu- 
tion decays as a (noninteger) power of geodesic distance 
to the singularity. For a recent review of critical collapse, 
see |. 

2. Beyond spherical symmetry 

An obvious question is if critical collapse is tied to 
spherical symmetry. Abrahams and Evans 0] fine-tuned 
the collapse of axisymmetric gravitational waves, and 
found evidence for black hole mass scaling, universality, 
and a self-similar critical solution. This shows that crit- 
ical phenomena can occur for axisymmetric, highly non- 
spherical initial data. In a complementary approach, we 
take here a known critical solution in spherical symme- 
try, and perturb it non-spherically. We find that the only 
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growing mode is the previously known spherical one. 

This result is weaker than that of Abrahams and Evans 
in being only perturbative. But a linearized result es- 
tablishes a result in an open neighborhood, up to lin- 
earization stability. In spherical symmetry the basin of 
attraction of the critical solution was empirically found 
to be the entire black hole threshold, even far from the 
critical solution This makes nonlinear instabilities 
seem extremely unlikely, and suggests that the open set 
may be in fact quite large. Generally it is remarkable 
how much in critical collapse can be understood quan- 
titatively through linear perturbation calculations: not 
only the critical exponent 7, but also universality with 
respect to initial data and matter models, and critical 
exponents for the black hole electric charge and angular 
momentum. 

Our result is stronger than that of Abrahams and 
Evans in not being limited to axisymmetry. Choptuik's 
spherically symmetric results showed that a genericity 
condition on initial data is necessary in formulating cos- 
mic censorship. Here we have demonstrated for the 
first time that this condition is sharp: The set of "non- 
generic" initial data that form a naked singularity can 
locally have codimension one in the full phase space. 



3. The background solution 

We write the general spherically symmetric spacetimc 
as a manifold M = M 2 x S* 2 with metric 



g av = diag (g AB 



r 2 lab) 



(1) 



where gAB is an arbitrary metric on M 2 , r is a scalar 
on M 2 , with r = defining the boundary of M 2 , and 
7ab is the unit curvature metric on S 2 . This spacetime is 
continously self-similar if there are coordinates x and r 
on M 2 such that in these coordinates the rescaled met- 
ric coefficients e 2T gAB and e 2r r 2 do not depend on r. 
It is discretely self-similar if they are periodic in r. A 
linearized spherical perturbation of a continously self- 
similar spacetime can be decomposed into modes of the 
schematic form e Xr f(x). The critical solution has pre- 
cisely one mode with ReX > 0. This is the mode that has 
to be beaten down by fine-tuning, while all other modes 
die away naturally as the curvature singularity t = 00 is 
approached. 

As a matter model we choose the perfect fluid with 
equation of state p = c 2 p, with c 2 a constant. While our 
equations hold for < c 2 < 1, we have carried out the 
numerical calculations only for c 2 = 1/3, the equation of 
state of a radiation gas. 

The coordinates x and r adapted to self-similarity are 
not unique, and they need not be fixed for the purpose of 
the following discussion. In our numerical calculations, 
however, we make a coordinate choice that is based on the 
Schwarzschild-like form ds 2 = -N 2 dt 2 + A 2 dr 2 + r 2 dtl 2 



of the metric. The coordinate transformation t = —e 
r = sxe~ T brings this into the form 



2t ( -N 2 + s 2 x 2 A 2 ~s 2 xA 2 
Symm 



9AB- e \ S 2 A 2 



2 -2t 2 2 

r = e s x . 



(2) 
(3) 



As a final gauge condition, we impose N = 1 at x = 
0, that is, t is central proper time. The spacetime in 
these coordinates is self-similar if N(x,t) and A(x, r) 
are functions of x only. 

The background stress-energy tensor is 



where 



i M „ = diag (t A B,c 2 s pr 2 "j ab ) 



tAB = (1 + c 2 )pu A u B + c 2 s pg A B, 



(4) 



(5) 



and — (ua , 0) is the fluid 4- velocity. The calcula- 
tion of the background solution as a non-linear boundary 
value problem was carried out along the lines of 
For the density we make the ansatz p — e 2T p(x), and 
for the radial velocity, u a t.a = v(x). Under this ansatz 
the Einstein and matter equations go over into a system 
of coupled ordinary differential equations (ODEs). These 
have a regular singular point (sonic point) where the sur- 
face x — const becomes an ingoing matter characteristic. 
The solution we want is the one uniquely specified by 
regularity at the origin x — and at the sonic point. 
We define the constant s so that the sonic point is at 
x = 1 , and solve the ODEs as a boundary value problem 
between x = and x = 1, solving for s as a nonlinear 
eigenvalue. 



4- Perturbation method 



Now one could perturb the critical solution Z(x) with 



an ansatz AZ(i,t) 



f(x) and solve the linearized 



boundary value problem for the discrete (complex) eigen- 
values A and mode functions f{x). This program has 
been carried out both for continously and discretely 
H self-similar solutions, and has allowed a precise numer- 
ical calculation of the critical exponent 7, which can be 
shown by dimensional analysis to be simply the inverse of 
the one positive real eigenvalue A. As a byproduct, one 
can determine the entire perturbation spectrum (which 
is discrete). 

Here, we mainly want to know is if there is any eigen- 
value with positive real part among the non-spherical 
perturbations, besides the one in the spherical pertur- 
bations. To answer this yes/no question, we write the 
perturbation equations as evolution equations in the time 
variable r. We decompose into spherical harmonics, and 
consider each value of I and m separately. Because the 
background is spherically symmetric, the dynamics of 
perturbations are the same for all values of m (for given 
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I). We evolve generic initial data for these equations for 
a sufficiently large interval of r. Generic data, with no 
field vanishing anywhere, constitute a superposition of 
all the (unknown) perturbation modes. 

In the time evolution, the mode with the largest ReX 
takes over after a transition period, and both that A and 
the corresponding f(x) can be simply read off from the 
late-time data 0. For I = in particular, this allows us 
to check our procedure by reading off the critical expo- 
nent as 7 ~ 0.36 in agreement with previous calculations 

§0- 

On a numerical grid, the frequency of modes that can 
be represented is limited by the grid spacing, so that 
we are only probing a large finite-dimensional subspace 
of all possible modes. Nevertheless, one can rule out 
the existence of unstable modes at very high spatial fre- 
quency (with respect to the coordinate x) by the fol- 
lowing argument. The perturbation equations form a 
system of linear wave equations with x-dependent coeffi- 
cients and an x-dependent mass matrix. High-frequency 
modes propagate essentially by geometric optics, and the 
mass terms are irrelevant for them. Therefore all high- 
frequency modes have the same dynamics. If they are 
shown to be decaying at frequencies still resolved by the 
numerical grid but at which the mass-like terms can al- 
ready be neglected with respect to derivative terms, then 
we can be sure that even higher frequencies not resolved 
on the numerical grid will decay in the same way. The 
geometric optics argument also applies to large values of 
I, which labels spatial frequencies in the angular coor- 
dinates 9 and tp. However, high I modes actually decay 
faster in r with increasing I, by a factor of e~ lT . Roughly 
speaking, this is because regular perturbations must be 
of 0(r l ) at the origin. 

The spherical perturbation equations are obtained sim- 
ply by linearizing the spherical field equations. There 
is no difficulty in writing them as evolution equations 
in time coordinate r and radial coordinate x. Bound- 
ary conditions at x = arise from demanding regularity 
at the center of spherical symmetry (all fields must be 
either even or odd in x). The surface x — 1 is an ingo- 
ing characteristic of the spherical perturbation equations, 
that is, a sound cone. No physical boundary conditions 
for perturbations are required there. Numerically, this 
absence of boundary conditions is implemented by a fi- 
nite difference scheme that is aware of the characteristics. 
Then we can bound the numerical domain by x — on 
one side and x = 1, or any larger constant x, on the 
other and evolve to arbitrary values of t. If one allows 
for non-spherical perturbations, one encounters gravita- 
tional as well as sound waves. The numerical domain 
for the entire system of perturbations then has to be ex- 
tended to the ingoing light cone, that is to x = x c defined 
by sx c A(x c ) = N(x c ). The background solution is eas- 
ily extended from x — 1 to x — x c and beyond. Fig. 
1 illustrates the coordinates, the characteristics, and the 
numerical domain. 



5. Gauge-invariant perturbations 

In going beyond spherical symmetry, one also has to 
deal with gauge-dependence and the presence of con- 
straints in the linearized Einstein and matter equations. 
Throughout this letter, we use the formalism and no- 
tation of Gerlach and Sengupta (GS) |fj"ofl . Any linear 
perturbation around spherical symmetry can be decom- 
posed into scalar, vector and tensor fields on M 2 times 
spherical harmonics Y™ on S 2 . Different l,m decouple. 
In the following we consider one value of (,ra at a time, 
and no longer write these indices on Y. Spherical har- 
monic vector fields on S 2 are Y a and S a = £a b Y. b , where a 
colon indicates the covariant derivative on S 2 , j a b-.c = 0, 
and t a b is the covariantly constant totally antisymmet- 
ric tensor on S 2 , e a b-.c = 0. Tensor perturbations based 
on Y and its derivatives (even or polar perturbations) 
decouple from tensor perturbations based on S a and its 
derivatives (odd or axial perturbations) . In the following 
we consider even and odd perturbations separately. 

GS express the 10 metric perturbations in the 2 + 2 
split. They calculate their transformation to linear order 
under the 4 infinitesimal coordinate changes, and find 6 
linear combinations of metric perturbations (with coef- 
ficients depending on the background) that are gauge- 
invariant to this order. The 10 stress-energy perturba- 
tions can be combined to form 10 gauge- invariant ones. 
The linearized Einstein equations can then be expressed 
in terms of those 16 variables alone. The remaining 4 
variables are pure gauge in the sense that one can give 
them arbitrary values, and reconstruct all 20 metric and 
stress-energy perturbations in the particular gauge de- 
termined that way. 

Calculation with, and interpretation of, the GS vari- 
ables is simplified by the fact that there is a gauge in 
which they can be identified directly with 16 gauge- 
dependent perturbations, while the remaining 4 gauge- 
dependent metric perturbations vanish. This gauge is the 
Regge- Wheeler gauge. In order to keep the notation sim- 
ple, we present the 16 gauge-invariant perturbations in 
this manner. The split between the odd and even sectors 
is as follows. There are 3 odd metric perturbations, and 
1 odd infinitesimal coordinate transformation, leaving 2 
odd gauge-invariant perturbations, in the form of a vec- 
tor field on M 2 . There are 7 — 3 = 4 even gauge-invariant 
metric perturbations, in the form of a symmetric tensor 
and a scalar. There are 3 odd and 7 even gauge-invariant 
matter perturbations. The general odd metric and mat- 
ter perturbations, in Regge- Wheeler gauge but expressed 
through the gauge-invariants, are 

^^Symm J' (6) 

At -( ° LaY " ) (7) 

^Symm L(S a:b + S b:a ) J' {t) 

The general even metric and matter perturbations are 
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At - ( Tab TaY:1 > \ (Q\ 

\ Symm TV 2 7afc + T 2 Y ab J ' W 

We now define gauge-invariant fluid perturbations, and 
present them once more by giving the equivalent gauge- 
dependent perturbations in Regge- Wheeler gauge. The 
odd perturbation of the fluid 4-velocity is 

Au„ = (0,pS a ). (10) 

(3 parameterizes axial fluid rotation. The even perturba- 
tion of the fluid 4-velocity is 

Au M = {Au A , aY. a ), where Au A = jn A + k AB u B . (11) 

Here n A = e AB ug, with c A b the totally antisymmetric 
covariant unit tensor on M 2 . Note that u A Au A = 0. 
a parameterizes axial fluid motion, while 7 parameter- 
izes perturbations of the radial fluid motion. The den- 
sity perturbation also belongs to the even sector, and 
we parameterize it as Ap = luY p. (By virtue of our 
equation of state, Ap = c? s Ap.) The general gauge- 
invariant stress-energy perturbations of GS are related 
to the gauge-invariant fluid perturbations as follows. In 
the odd sector, we have 

L A = (3{l + c 2 s )pu A , L = 0, (12) 

and in the even sector, 

T A = a(l + c 2 )pu A , T 1 = (lo + k)c 2 p, T 2 = 0, (13) 

Tab = Lut A B + 2Au( A u B )(l + c? s )p + kABC 2 s p. (14) 



6. Odd perturbations 

We must now extract a well-posed initial value problem 
from the gauge- invariant perturbation equations. For the 
odd sector this is straightforward. The one nontrivial 
matter conservation equation is {t 2 La)^ a — (where 
9ab\c = 0), or 

{fir 2 pu% A = 0. (15) 

This equation is an advection equation for 0, and can be 
solved independently of all other perturbations, from 
given on an initial spacelike hypersurface of M 2 . GS have 
shown that by defining the scalar II = t AE ' {r~ 2 k A )\B the 
Einstein equations for k A can always be reduced to the 
scalar wave equation 

[r- 2 (r 4 U) lA ]\ A -(1 + 2) (J - 1)11 = 16ne AB L AlBl (16) 

which they call the odd-parity master equation. The full 
k A can be reconstructed by quadratures once this equa- 
tion has been solved for II. (Note that the source L A is 
already known in the case of perfect fluid matter.) 



7. Even perturbations 

The even perturbations are more entangled. For a we 
once again have an advection equation, but now with 
sources. The first-order equations for the density per- 
turbation lo and radial velocity perturbation 7 can be 
combined to form a single wave equation at the speed 
of sound for either uj or 7. For fixed metric perturbation 
k A B and k, the initial value problem would then be clear. 
Unfortunately, one apparently cannot extract a master 
equation for the metric perturbations from the linearized 
Einstein equations for arbitrary matter, but one always 
has a system with constraints. Instead, we follow Seidel 
p"l| in first focussing attention on those components of 
the linear Einstein equations with vanishing matter sides. 
4 such components exist, because the 7 even stress-energy 
perturbations are linear in only the 3 even matter per- 
turbations. One of these is 

k/ = -I671-T 2 = 0. (17) 

It is natural to decompose the trace-free tensor k A B co- 
variantly into two scalars, via 

k A B = 4>{u A u B + n A n B ) + ^{u A n B + n A u B ), (18) 

and to introduce the frame derivatives 

/ = u A f lA , f = n A f [A . (19) 

The three remaining source-free Einstein equations can 
then be written as 

-(Xy+(X')'=S X , (20) 

-(ky+c 2 s (k'y=s k , (21) 

-4= V (22) 

where \ = 4> — k replaces <f>. The source terms S x , Sk 
and S^p are linear in x, fc, ip and their first derivatives 
x', k', ip', x, and k, but do not contain -0. While the 
highest derivatives of \ form a wave equation with char- 
acteristics given by the metric g A s = —U A Ub + n A ns, k 
obeys a wave equation with characteristics given by the 
"fluid metric" — u a ub + c~ 2 n A riB- These characteristics 
have speed c s relative to the fluid. Finally, ip is advected 
with the fluid. Therefore x characterizes gravitational 
waves, k sound waves, and ip polar fluid flow. The met- 
ric perturbation k "knows about" the speed of sound be- 
cause we have used Ap = dp/dpAp = c 2 Ap in finding 
the source-free linearized Einstein equations. If we add 
to these three equations the identity {f'u A )\ A — (fn A )\ A 
and the definitions (ff9|), we have a complete first-order 
system of equations. The variables %, x, k, k and ip can 
be set freely on a spacelike hypersurface in M 2 . From the 
equations of motion one can see that any regular solution 
must scale at r = as k ~ r , ip ~ r l+1 , and x ~ r l+2 . 
This follows also from the requirement that the metric 
perturbation (^) be a regular tensor in four dimensions 
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at r = 0. Furthermore, the perturbed metric remains 
continuously self-similar if k, ip and x are independent of 
r. The perturbed metric is discretely self-similar if these 
fields are periodic in r. 

The Einstein equations we have not used yet give the 
matter perturbations directly in terms of the metric per- 
turbations. As a check of the correctness of our equations 
and their numerical implementation, we have numerically 
differentiated the numerical solution and verified that the 
perturbed matter equations of motion (or Bianchi iden- 
tities) are obeyed. 

The case of I = 1 even perturbations is not covered by 
the framework of GS and has to be treated separately. 
Clearly / = 1 does not admit gravitational waves, so that 
we can use Newtonian intuition. The I — 1 matter per- 
turbations are pure gauge, corresponding to an initial 
displacement and velocity of the spherical background 
solution. 



8. Numerical method 

We conclude with a remark on the numerical imple- 
mentation. The linearized field equations are of the form 



Bu. 



(23) 



du ^du 
dr dx 

In order to make the numerical evolution stable even 
though the lines of constant x change from timelike to 
spacelike in our computational domain, it is essential to 
use a characteristic scheme jl^] . Let V be the matrix of 
(column) eigenvectors of A. Let A be the diagonal ma- 
trix composed of the corresponding eigenvalues. Then 
A = VAV^ 1 . Let A+ be the diagonal matrix with zeros 
in the place of the negative eigenvalues. Define A_, A + 
and A- in the obvious manner, so that A = A + + A_. 
We have used the numerical scheme 



At 



A 4 



Aa; 



A^ 



*3-l 



Aa 



Bu-, 



(24) 

which is first-order accurate, and stable even for super- 
luminal shift. The matrix A is just sparse enough for its 
eigenvalues and eigenvectors to be calculated in closed 
form. As expected, the characteristics are the fluid world 
lines, light cones and sound cones. 
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FIG. 1. The horizontal thin lines accumulating at t = are 
lines of constant r. The fanning thin lines are lines of constant 
x. The inner dashed line is the ingoing sound cone. The 
outer dashed line is the ingoing light cone. The dark shaded 
region is the numerical domain of the background, spherical 
and odd perturbations calculations. The total shaded region 
is the numerical domain of the even perturbations calculation. 
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